#################################################################################
###Replication Files for Coalition Quality and Multinational Conflict Outcomes###
###By Skyler J. Cranmer and Elizabeth J. Menninga                             ###
###For publication at International Interactions                              ###
###Last Updated 7/12/2017                                                     ###
#################################################################################

library(MCMCpack)
library(xtable)
library(ROCR)
library(PRROC)
library(separationplot)

set.seed(235813)

###Functions for formatting table output & used to extract quantities of interest
interpret <- function(mat){
  prob.effect <- function(x){
    x[x > 0] <- 1
    x[x < 0] <- 0
    out <- c(mean(x), 1-mean(x))
    return(out)
  }
  out <- matrix(NA, nrow=ncol(mat), ncol=4)
  for(i in 1:ncol(mat)){
    out[i,1] <- mean(mat[,i])
    out[i,2] <- sd(mat[,i])
    out[i,3] <- prob.effect(mat[,i])[1]
    out[i,4] <- prob.effect(mat[,i])[2]
  }
  rownames(out) <- colnames(mat)
  colnames(out) <- c("Mean","SD","Prob > 0","Prob < 0")
  out <- round(out,2)
  return(out)
}

tex.interpret <- function(mat){
  prob.effect <- function(x){
    x[x > 0] <- 1
    x[x < 0] <- 0
    out <- c(mean(x), 1-mean(x))
    return(out)
  }
  out <- matrix(NA, nrow=ncol(mat), ncol=4)
  for(i in 1:ncol(mat)){
    out[i,1] <- mean(mat[,i])
    out[i,2] <- sd(mat[,i])
    out[i,3] <- prob.effect(mat[,i])[1]
    out[i,4] <- prob.effect(mat[,i])[2]
  }
  rownames(out) <- colnames(mat)
  colnames(out) <- c("Mean","SD","Prob > 0","Prob < 0")
  out <- round(out,2)
  col1 <- paste(out[,1]," (",out[,2],") ", sep="" )
  out2 <- cbind(col1, out[,3:4])
  colnames(out2) <- c("Mean (SD)", "Prob > 0","Prob < 0")
  return(out2)
}

###Table 1

#Set working directory
coal<-read.csv("Replication Files/coal.csv")

##Variables are in order as they are reported in the table
##ap=proportion experienced; awin=success rate; aally=proportion allied; aq=previous collaboration;
##qsuccess=collaborative success rate; ademvs=political variance; lnuma=log size Side A; caprat2=capability ratio;
##lnumb=log size Side B; adistm=mean geographic distance
formcoal <- as.formula("aout_vic ~ ap + awin  + aally + aq + ademvs + lnuma + caprat2 + lnumb+ adistm")
coal_est <- MCMClogit(formcoal, burnin=1000000, mcmc=10000000, data=coal)

formcoal2 <- as.formula("aout_vic ~ ap  + aally + aq + qsuccess +  ademvs + lnuma + caprat2 + lnumb+ adistm")
coal_est2 <- MCMClogit(formcoal2, burnin=1000000, mcmc=10000000, data=coal)

xtable(cbind(tex.interpret(coal_est),tex.interpret(coal_est2)))

##Table 2
av<-read.csv("Replication Files/av.csv") #Initiating coalitions data
dv<-read.csv("Replication Files/dv.csv") #Defending coalitions data
av$caprat2 <- av$acap / (av$acap + av$bcap)
dv$caprat2 <- dv$acap / (dv$acap + dv$bcap)

#Initiating models; variables are same as above
formBaall1 <- as.formula("av_out_all ~ ap + awin + aally + aq +  ademvs+ l.numa + caprat2 + l.numb+ adistm")
formBaall2<- as.formula("av_out_all ~ ap +  aally + aq + qsuccess + ademvs+ l.numa + caprat2 + l.numb+ adistm")
asab_all <- MCMClogit(formBaall1, burnin=1000000, mcmc=5000000, data=av)
asab_all2 <- MCMClogit(formBaall2, burnin=1000000, mcmc=5000000, data=av)

#Defending models
formBdall1 <- as.formula("dv_out_all ~ bp + bwin + bally + bq +  bdemvs+ l.numa + caprat2 + l.numb+ bdistm")
formBdall2 <- as.formula("dv_out_all ~ bp +  bally + bq + qsuccess + bdemvs+ l.numa + caprat2 + l.numb+ bdistm")
asdb_all <- MCMClogit(formBdall1, burnin=1000000, mcmc=5000000, data=dv)
asdb_all2 <- MCMClogit(formBdall2, burnin=1000000, mcmc=5000000, data=dv)

xtable(cbind(tex.interpret(asab_all),tex.interpret(asab_all2), tex.interpret(asdb_all), tex.interpret(asdb_all2)))

########################## Figure 1
pprobs <- function(posterior, X){
  Xb <- as.matrix(X) %*% t(posterior)
  probs <- 1/(1+exp(-Xb))
  return(probs)
}

# Taking a sample from the posterior distributions 
All_sam<-coal_est[sample(nrow(coal_est),size=10000,replace=F),] 

# ap
# Create Data frame that captures variation in ap with all others at their means #
ap_range<-seq(0, max(coal$ap,na.rm=T)*1.1, length.out=20)
ap_pp_plot<-data.frame(rep(1, length.out=20), ap_range, rep(mean(coal$awin), length.out=20), rep(mean(coal$aally), length.out=20),rep(mean(coal$aq), length.out=20),rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
ap_pp_plot_pred<-pprobs(ap_pp_plot, All_sam)
ap_mean<-apply(ap_pp_plot_pred, 2, mean)
ap_ci<-apply(ap_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# awin
# Create Data frame that captures variation in awin with all others at their means #
awin_range<-seq(0, max(coal$awin,na.rm=T)*1.1, length.out=20)
awin_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), awin_range,  rep(mean(coal$aally), length.out=20),rep(mean(coal$aq), length.out=20),rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
awin_pp_plot_pred<-pprobs(awin_pp_plot, All_sam)
awin_mean<-apply(awin_pp_plot_pred, 2, mean)
awin_ci<-apply(awin_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# aally
# Create Data frame that captures variation in aally with all others at their means #
aally_range<-seq(0, max(coal$aally,na.rm=T)*1.1, length.out=20)
aally_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$awin), length.out=20), aally_range,  rep(mean(coal$aq), length.out=20),rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
aally_pp_plot_pred<-pprobs(aally_pp_plot, All_sam)
aally_mean<-apply(aally_pp_plot_pred, 2, mean)
aally_ci<-apply(aally_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# aq
# Create Data frame that captures variation in aq with all others at their means #
aq_range<-seq(0, max(coal$aq,na.rm=T)*1.1, length.out=20)
aq_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$awin), length.out=20), rep(mean(coal$aally), length.out=20), aq_range, rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
aq_pp_plot_pred<-pprobs(aq_pp_plot, All_sam)
aq_mean<-apply(aq_pp_plot_pred, 2, mean)
aq_ci<-apply(aq_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# ademvs
# Create Data frame that captures variation in ademvs with all others at their means #
ademvs_range<-seq(0, max(coal$ademvs,na.rm=T)*1.1, length.out=20)
ademvs_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$awin), length.out=20), rep(mean(coal$aally), length.out=20), rep(mean(coal$aq), length.out=20), ademvs_range, rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
ademvs_pp_plot_pred<-pprobs(ademvs_pp_plot, All_sam)
ademvs_mean<-apply(ademvs_pp_plot_pred, 2, mean)
ademvs_ci<-apply(ademvs_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# lnuma
# Create Data frame that captures variation in lnuma with all others at their means #
lnuma_range<-seq(0, max(coal$lnuma,na.rm=T)*1.1, length.out=20)
lnuma_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$awin), length.out=20), rep(mean(coal$aally), length.out=20), rep(mean(coal$aq), length.out=20), rep(mean(coal$ademvs), length.out=20), lnuma_range, rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
lnuma_pp_plot_pred<-pprobs(lnuma_pp_plot, All_sam)
lnuma_mean<-apply(lnuma_pp_plot_pred, 2, mean)
lnuma_ci<-apply(lnuma_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

windows()
par(mfrow=c(2,3)) 

plot(c(0,max(coal$ap,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Experienced") 
lines(ap_range, ap_mean, lwd=3, lty=1)
lines(ap_range, ap_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(ap_range, ap_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(ap_range, ap_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(ap_range, ap_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$awin,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Success Rate") 
lines(awin_range, awin_mean, lwd=3, lty=1)
lines(awin_range, awin_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(awin_range, awin_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(awin_range, awin_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(awin_range, awin_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$aally,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Allied") 
lines(aally_range, aally_mean, lwd=3, lty=1)
lines(aally_range, aally_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(aally_range, aally_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(aally_range, aally_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(aally_range, aally_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$aq,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Previous Collaboration") 
lines(aq_range, aq_mean, lwd=3, lty=1)
lines(aq_range, aq_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(aq_range, aq_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(aq_range, aq_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(aq_range, aq_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$ademvs,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Political Variance") 
lines(ademvs_range, ademvs_mean, lwd=3, lty=1)
lines(ademvs_range, ademvs_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(ademvs_range, ademvs_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(ademvs_range, ademvs_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(ademvs_range, ademvs_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$lnuma,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Log Coalition Size") 
lines(lnuma_range, lnuma_mean, lwd=3, lty=1)
lines(lnuma_range, lnuma_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(lnuma_range, lnuma_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(lnuma_range, lnuma_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(lnuma_range, lnuma_ci[3,],  lwd=1.5, lty=1, col="grey50")

######################## Figure 2

# Taking a sample from the posterior distributions 
All_sam<-coal_est2[sample(nrow(coal_est2),size=10000,replace=F),] 

# ap
# Create Data frame that captures variation in ap with all others at their means #
ap_range<-seq(0, max(coal$ap,na.rm=T)*1.1, length.out=20)
ap_pp_plot<-data.frame(rep(1, length.out=20), ap_range, rep(mean(coal$aally), length.out=20),rep(mean(coal$aq), length.out=20), rep(mean(coal$qsuccess), length.out=20), rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
ap_pp_plot_pred<-pprobs(ap_pp_plot, All_sam)
ap_mean<-apply(ap_pp_plot_pred, 2, mean)
ap_ci<-apply(ap_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# aally
# Create Data frame that captures variation in aally with all others at their means #
aally_range<-seq(0, max(coal$aally,na.rm=T)*1.1, length.out=20)
aally_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), aally_range,  rep(mean(coal$aq), length.out=20), rep(mean(coal$qsuccess), length.out=20), rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
aally_pp_plot_pred<-pprobs(aally_pp_plot, All_sam)
aally_mean<-apply(aally_pp_plot_pred, 2, mean)
aally_ci<-apply(aally_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# aq
# Create Data frame that captures variation in aq with all others at their means #
aq_range<-seq(0, max(coal$aq,na.rm=T)*1.1, length.out=20)
aq_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$aally), length.out=20), aq_range, rep(mean(coal$qsuccess), length.out=20), rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
aq_pp_plot_pred<-pprobs(aq_pp_plot, All_sam)
aq_mean<-apply(aq_pp_plot_pred, 2, mean)
aq_ci<-apply(aq_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# qsuccess
# Create Data frame that captures variation in qsuccess with all others at their means #
qsuccess_range<-seq(0, max(coal$qsuccess,na.rm=T)*1.1, length.out=20)
qsuccess_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$aally), length.out=20),rep(mean(coal$aq), length.out=20), qsuccess_range,  rep(mean(coal$ademvs), length.out=20),rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
qsuccess_pp_plot_pred<-pprobs(qsuccess_pp_plot, All_sam)
qsuccess_mean<-apply(qsuccess_pp_plot_pred, 2, mean)
qsuccess_ci<-apply(qsuccess_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# ademvs
# Create Data frame that captures variation in ademvs with all others at their means #
ademvs_range<-seq(0, max(coal$ademvs,na.rm=T)*1.1, length.out=20)
ademvs_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$aally), length.out=20), rep(mean(coal$aq), length.out=20), rep(mean(coal$qsuccess), length.out=20), ademvs_range, rep(mean(coal$lnuma), length.out=20),rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
ademvs_pp_plot_pred<-pprobs(ademvs_pp_plot, All_sam)
ademvs_mean<-apply(ademvs_pp_plot_pred, 2, mean)
ademvs_ci<-apply(ademvs_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

# lnuma
# Create Data frame that captures variation in lnuma with all others at their means #
lnuma_range<-seq(0, max(coal$lnuma,na.rm=T)*1.1, length.out=20)
lnuma_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(coal$ap), length.out=20), rep(mean(coal$aally), length.out=20), rep(mean(coal$aq), length.out=20), rep(mean(coal$qsuccess), length.out=20), rep(mean(coal$ademvs), length.out=20), lnuma_range, rep(mean(coal$caprat2), length.out=20),rep(mean(coal$lnumb), length.out=20),rep(mean(coal$adistm), length.out=20))

# Create pred. probs. 
lnuma_pp_plot_pred<-pprobs(lnuma_pp_plot, All_sam)
lnuma_mean<-apply(lnuma_pp_plot_pred, 2, mean)
lnuma_ci<-apply(lnuma_pp_plot_pred, 2, quantile, probs=c(0.05, 0.16, 0.84, 0.95))

windows()
par(mfrow=c(2,3)) 

plot(c(0,max(coal$ap,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Experienced") 
lines(ap_range, ap_mean, lwd=3, lty=1)
lines(ap_range, ap_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(ap_range, ap_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(ap_range, ap_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(ap_range, ap_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$aally,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Allied") 
lines(aally_range, aally_mean, lwd=3, lty=1)
lines(aally_range, aally_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(aally_range, aally_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(aally_range, aally_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(aally_range, aally_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$aq,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Previous Collaboration") 
lines(aq_range, aq_mean, lwd=3, lty=1)
lines(aq_range, aq_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(aq_range, aq_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(aq_range, aq_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(aq_range, aq_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$qsuccess,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Collaborative Success Rate") 
lines(qsuccess_range, qsuccess_mean, lwd=3, lty=1)
lines(qsuccess_range, qsuccess_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(qsuccess_range, qsuccess_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(qsuccess_range, qsuccess_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(qsuccess_range, qsuccess_ci[3,],  lwd=1.5, lty=1, col="grey50")

plot(c(0,max(coal$ademvs,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Political Variance") 
lines(ademvs_range, ademvs_mean, lwd=3, lty=1)
lines(ademvs_range, ademvs_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(ademvs_range, ademvs_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(ademvs_range, ademvs_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(ademvs_range, ademvs_ci[3,],  lwd=1.5, lty=1, col="grey50")


plot(c(0,max(coal$lnuma,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Log Coalition Size") 
lines(lnuma_range, lnuma_mean, lwd=3, lty=1)
lines(lnuma_range, lnuma_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(lnuma_range, lnuma_ci[4,],  lwd=1.5, lty=2, col="grey50")
lines(lnuma_range, lnuma_ci[2,], lwd=1.5, lty=1, col="grey50")
lines(lnuma_range, lnuma_ci[3,],  lwd=1.5, lty=1, col="grey50")


############################## Figure 3

# Taking a sample from the posterior distributions for initiating coalitions
All_sam<-asab_all[sample(nrow(asab_all),size=10000,replace=F),] 

# ap
# Create Data frame that captures variation in ap with all others at their means #
ap_range<-seq(0, max(av$ap,na.rm=T)*1.1, length.out=20)
ap_pp_plot<-data.frame(rep(1, length.out=20), ap_range, rep(mean(av$awin), length.out=20), rep(mean(av$aally), length.out=20),rep(mean(av$aq), length.out=20), rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
ap_pp_plot_pred<-pprobs(ap_pp_plot, All_sam)
ap_mean<-apply(ap_pp_plot_pred, 2, mean)
ap_ci<-apply(ap_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# awin
# Create Data frame that captures variation in awin with all others at their means #
awin_range<-seq(0, max(av$awin,na.rm=T)*1.1, length.out=20)
awin_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), awin_range, rep(mean(av$aally), length.out=20),rep(mean(av$aq), length.out=20), rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
awin_pp_plot_pred<-pprobs(awin_pp_plot, All_sam)
awin_mean<-apply(awin_pp_plot_pred, 2, mean)
awin_ci<-apply(awin_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# aally
# Create Data frame that captures variation in aally with all others at their means #
aally_range<-seq(0, max(av$aally,na.rm=T)*1.1, length.out=20)
aally_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$awin), length.out=20), aally_range,  rep(mean(av$aq), length.out=20), rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
aally_pp_plot_pred<-pprobs(aally_pp_plot, All_sam)
aally_mean<-apply(aally_pp_plot_pred, 2, mean)
aally_ci<-apply(aally_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# aq
# Create Data frame that captures variation in aq with all others at their means #
aq_range<-seq(0, max(av$aq,na.rm=T)*1.1, length.out=20)
aq_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$awin), length.out=20), rep(mean(av$aally), length.out=20), aq_range, rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
aq_pp_plot_pred<-pprobs(aq_pp_plot, All_sam)
aq_mean<-apply(aq_pp_plot_pred, 2, mean)
aq_ci<-apply(aq_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# ademvs
# Create Data frame that captures variation in ademvs with all others at their means #
ademvs_range<-seq(0, max(av$ademvs,na.rm=T)*1.1, length.out=20)
ademvs_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$awin), length.out=20), rep(mean(av$aally), length.out=20), rep(mean(av$aq), length.out=20), ademvs_range, rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
ademvs_pp_plot_pred<-pprobs(ademvs_pp_plot, All_sam)
ademvs_mean<-apply(ademvs_pp_plot_pred, 2, mean)
ademvs_ci<-apply(ademvs_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# l.numa
# Create Data frame that captures variation in l.numa with all others at their means #
lnuma_range<-seq(0, max(av$l.numa,na.rm=T)*1.1, length.out=20)
lnuma_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$awin), length.out=20), rep(mean(av$aally), length.out=20), rep(mean(av$aq), length.out=20), rep(mean(av$ademvs), length.out=20), lnuma_range, rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
lnuma_pp_plot_pred<-pprobs(lnuma_pp_plot, All_sam)
lnuma_mean<-apply(lnuma_pp_plot_pred, 2, mean)
lnuma_ci<-apply(lnuma_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))


# Taking a sample from the posterior distributions for defending coalitions
All_sam_def<-asdb_all[sample(nrow(asdb_all),size=10000,replace=F),] 

# bp
# Create Data frame that captures variation in bp with all others at their means #
bp_range<-seq(0, max(dv$bp,na.rm=T)*1.1, length.out=20)
bp_pp_plot<-data.frame(rep(1, length.out=20), bp_range, rep(mean(dv$bwin), length.out=20), rep(mean(dv$bally), length.out=20),rep(mean(dv$bq), length.out=20), rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bp_pp_plot_pred<-pprobs(bp_pp_plot, All_sam_def)
bp_mean<-apply(bp_pp_plot_pred, 2, mean)
bp_ci<-apply(bp_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# bwin
# Create Data frame that captures variation in bwin with all others at their means #
bwin_range<-seq(0, max(dv$bwin,na.rm=T)*1.1, length.out=20)
bwin_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), bwin_range,  rep(mean(dv$bally), length.out=20), rep(mean(dv$bq), length.out=20), rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bwin_pp_plot_pred<-pprobs(bwin_pp_plot, All_sam_def)
bwin_mean<-apply(bwin_pp_plot_pred, 2, mean)
bwin_ci<-apply(bwin_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# bally
# Create Data frame that captures variation in bally with all others at their means #
bally_range<-seq(0, max(dv$bally,na.rm=T)*1.1, length.out=20)
bally_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bwin), length.out=20), bally_range,  rep(mean(dv$bq), length.out=20), rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bally_pp_plot_pred<-pprobs(bally_pp_plot, All_sam_def)
bally_mean<-apply(bally_pp_plot_pred, 2, mean)
bally_ci<-apply(bally_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# bq
# Create Data frame that captures variation in bq with all others at their means #
bq_range<-seq(0, max(dv$bq,na.rm=T)*1.1, length.out=20)
bq_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bwin), length.out=20), rep(mean(dv$bally), length.out=20), bq_range, rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bq_pp_plot_pred<-pprobs(bq_pp_plot, All_sam_def)
bq_mean<-apply(bq_pp_plot_pred, 2, mean)
bq_ci<-apply(bq_pp_plot_pred, 2, quantile, probs=c(0.05,  0.95))

# bdemvs
# Create Data frame that captures variation in bdemvs with all others at their means #
bdemvs_range<-seq(0, max(dv$bdemvs,na.rm=T)*1.1, length.out=20)
bdemvs_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bwin), length.out=20), rep(mean(dv$bally), length.out=20), rep(mean(dv$bq), length.out=20), bdemvs_range, rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bdemvs_pp_plot_pred<-pprobs(bdemvs_pp_plot, All_sam_def)
bdemvs_mean<-apply(bdemvs_pp_plot_pred, 2, mean)
bdemvs_ci<-apply(bdemvs_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# l.numb
# Create Data frame that captures variation in l.numb with all others at their means #
lnumb_range<-seq(0, max(dv$l.numb,na.rm=T)*1.1, length.out=20)
lnumb_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bwin), length.out=20), rep(mean(dv$bally), length.out=20), rep(mean(dv$bq), length.out=20), rep(mean(dv$bdemvs), length.out=20), rep(mean(dv$l.numa), length.out=20), rep(mean(dv$caprat2), length.out=20), lnumb_range, rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
lnumb_pp_plot_pred<-pprobs(lnumb_pp_plot, All_sam_def)
lnumb_mean<-apply(lnumb_pp_plot_pred, 2, mean)
lnumb_ci<-apply(lnumb_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))


windows()
par(mfrow=c(2,3)) 

plot(c(0,max(av$ap,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Experienced") 
lines(ap_range, ap_mean, lwd=3, lty=1)
lines(ap_range, ap_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(ap_range, ap_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bp_range, bp_mean, lwd=3, lty=2)
lines(bp_range, bp_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bp_range, bp_ci[2,],  lwd=1.5, lty=2, col="grey50")


plot(c(0,max(av$awin,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Success Rate") 
lines(awin_range, awin_mean, lwd=3, lty=1)
lines(awin_range, awin_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(awin_range, awin_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bwin_range, bwin_mean, lwd=3, lty=2)
lines(bwin_range, bwin_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bwin_range, bwin_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$aally,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Allied") 
lines(aally_range, aally_mean, lwd=3, lty=1)
lines(aally_range, aally_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(aally_range, aally_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bally_range, bally_mean, lwd=3, lty=2)
lines(bally_range, bally_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bally_range, bally_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$aq,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Previous Collaboration") 
lines(aq_range, aq_mean, lwd=3, lty=1)
lines(aq_range, aq_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(aq_range, aq_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bq_range, bq_mean, lwd=3, lty=2)
lines(bq_range, bq_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bq_range, bq_ci[2,],  lwd=1.5, lty=2, col="grey50")


plot(c(0,max(av$ademvs,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Political Variance") 
lines(ademvs_range, ademvs_mean, lwd=3, lty=1)
lines(ademvs_range, ademvs_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(ademvs_range, ademvs_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bdemvs_range, bdemvs_mean, lwd=3, lty=2)
lines(bdemvs_range, bdemvs_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bdemvs_range, bdemvs_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$l.numa,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Log Coalition Size") 
lines(lnuma_range, lnuma_mean, lwd=3, lty=1)
lines(lnuma_range, lnuma_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(lnuma_range, lnuma_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(lnumb_range, lnumb_mean, lwd=3, lty=2)
lines(lnumb_range, lnumb_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(lnumb_range, lnumb_ci[2,],  lwd=1.5, lty=2, col="grey50")


############################## Figure 4

# Taking a sample from the posterior distributions for initiating coalitions
All_sam<-asab_all2[sample(nrow(asab_all2),size=10000,replace=F),] 

# ap
# Create Data frame that captures variation in ap with all others at their means #
ap_range<-seq(0, max(av$ap,na.rm=T)*1.1, length.out=20)
ap_pp_plot<-data.frame(rep(1, length.out=20), ap_range, rep(mean(av$aally), length.out=20),rep(mean(av$aq), length.out=20), rep(mean(av$qsuccess), length.out=20), rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
ap_pp_plot_pred<-pprobs(ap_pp_plot, All_sam)
ap_mean<-apply(ap_pp_plot_pred, 2, mean)
ap_ci<-apply(ap_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# aally
# Create Data frame that captures variation in aally with all others at their means #
aally_range<-seq(0, max(av$aally,na.rm=T)*1.1, length.out=20)
aally_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), aally_range,  rep(mean(av$aq), length.out=20), rep(mean(av$qsuccess), length.out=20), rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
aally_pp_plot_pred<-pprobs(aally_pp_plot, All_sam)
aally_mean<-apply(aally_pp_plot_pred, 2, mean)
aally_ci<-apply(aally_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# aq
# Create Data frame that captures variation in aq with all others at their means #
aq_range<-seq(0, max(av$aq,na.rm=T)*1.1, length.out=20)
aq_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$aally), length.out=20), aq_range, rep(mean(av$qsuccess), length.out=20), rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
aq_pp_plot_pred<-pprobs(aq_pp_plot, All_sam)
aq_mean<-apply(aq_pp_plot_pred, 2, mean)
aq_ci<-apply(aq_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# qsuccess
# Create Data frame that captures variation in qsuccess with all others at their means #
qsuccess_range<-seq(0, max(av$qsuccess,na.rm=T)*1.1, length.out=20)
qsuccess_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$aally), length.out=20),rep(mean(av$aq), length.out=20), qsuccess_range,  rep(mean(av$ademvs), length.out=20),rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
qsuccess_pp_plot_pred<-pprobs(qsuccess_pp_plot, All_sam)
qsuccess_mean<-apply(qsuccess_pp_plot_pred, 2, mean)
qsuccess_ci<-apply(qsuccess_pp_plot_pred, 2, quantile, probs=c(0.05,0.95))

# ademvs
# Create Data frame that captures variation in ademvs with all others at their means #
ademvs_range<-seq(0, max(av$ademvs,na.rm=T)*1.1, length.out=20)
ademvs_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$aally), length.out=20), rep(mean(av$aq), length.out=20), rep(mean(av$qsuccess), length.out=20), ademvs_range, rep(mean(av$l.numa), length.out=20),rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
ademvs_pp_plot_pred<-pprobs(ademvs_pp_plot, All_sam)
ademvs_mean<-apply(ademvs_pp_plot_pred, 2, mean)
ademvs_ci<-apply(ademvs_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# l.numa
# Create Data frame that captures variation in l.numa with all others at their means #
lnuma_range<-seq(0, max(av$l.numa,na.rm=T)*1.1, length.out=20)
lnuma_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(av$ap), length.out=20), rep(mean(av$aally), length.out=20), rep(mean(av$aq), length.out=20), rep(mean(av$qsuccess), length.out=20), rep(mean(av$ademvs), length.out=20), lnuma_range, rep(mean(av$caprat2), length.out=20),rep(mean(av$l.numb), length.out=20),rep(mean(av$adistm), length.out=20))

# Create pred. probs. 
lnuma_pp_plot_pred<-pprobs(lnuma_pp_plot, All_sam)
lnuma_mean<-apply(lnuma_pp_plot_pred, 2, mean)
lnuma_ci<-apply(lnuma_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))


# Taking a sample from the posterior distributions for defending coalitions
All_sam_def<-asdb_all2[sample(nrow(asdb_all2),size=10000,replace=F),] 

# bp
# Create Data frame that captures variation in bp with all others at their means #
bp_range<-seq(0, max(dv$bp,na.rm=T)*1.1, length.out=20)
bp_pp_plot<-data.frame(rep(1, length.out=20), bp_range, rep(mean(dv$bally), length.out=20),rep(mean(dv$bq), length.out=20), rep(mean(dv$qsuccess), length.out=20), rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bp_pp_plot_pred<-pprobs(bp_pp_plot, All_sam_def)
bp_mean<-apply(bp_pp_plot_pred, 2, mean)
bp_ci<-apply(bp_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# bally
# Create Data frame that captures variation in bally with all others at their means #
bally_range<-seq(0, max(dv$bally,na.rm=T)*1.1, length.out=20)
bally_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), bally_range,  rep(mean(dv$bq), length.out=20), rep(mean(dv$qsuccess), length.out=20), rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bally_pp_plot_pred<-pprobs(bally_pp_plot, All_sam_def)
bally_mean<-apply(bally_pp_plot_pred, 2, mean)
bally_ci<-apply(bally_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# bq
# Create Data frame that captures variation in bq with all others at their means #
bq_range<-seq(0, max(dv$bq,na.rm=T)*1.1, length.out=20)
bq_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bally), length.out=20), bq_range, rep(mean(dv$qsuccess), length.out=20), rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bq_pp_plot_pred<-pprobs(bq_pp_plot, All_sam_def)
bq_mean<-apply(bq_pp_plot_pred, 2, mean)
bq_ci<-apply(bq_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# qsuccess
# Create Data frame that captures variation in qsuccess with all others at their means #
bqsuccess_range<-seq(0, max(dv$qsuccess,na.rm=T)*1.1, length.out=20)
bqsuccess_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bally), length.out=20),rep(mean(dv$bq), length.out=20), qsuccess_range,  rep(mean(dv$bdemvs), length.out=20),rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bqsuccess_pp_plot_pred<-pprobs(bqsuccess_pp_plot, All_sam_def)
bqsuccess_mean<-apply(bqsuccess_pp_plot_pred, 2, mean)
bqsuccess_ci<-apply(bqsuccess_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# bdemvs
# Create Data frame that captures variation in bdemvs with all others at their means #
bdemvs_range<-seq(0, max(dv$bdemvs,na.rm=T)*1.1, length.out=20)
bdemvs_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bally), length.out=20), rep(mean(dv$bq), length.out=20), rep(mean(dv$qsuccess), length.out=20), bdemvs_range, rep(mean(dv$l.numa), length.out=20),rep(mean(dv$caprat2), length.out=20),rep(mean(dv$l.numb), length.out=20),rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
bdemvs_pp_plot_pred<-pprobs(bdemvs_pp_plot, All_sam_def)
bdemvs_mean<-apply(bdemvs_pp_plot_pred, 2, mean)
bdemvs_ci<-apply(bdemvs_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

# l.numb
# Create Data frame that captures variation in l.numb with all others at their means #
lnumb_range<-seq(0, max(dv$l.numb,na.rm=T)*1.1, length.out=20)
lnumb_pp_plot<-data.frame(rep(1, length.out=20), rep(mean(dv$bp), length.out=20), rep(mean(dv$bally), length.out=20), rep(mean(dv$bq), length.out=20), rep(mean(dv$qsuccess), length.out=20), rep(mean(dv$bdemvs), length.out=20), rep(mean(dv$l.numa), length.out=20), rep(mean(dv$caprat2), length.out=20), lnumb_range, rep(mean(dv$bdistm), length.out=20))

# Create pred. probs. 
lnumb_pp_plot_pred<-pprobs(lnumb_pp_plot, All_sam_def)
lnumb_mean<-apply(lnumb_pp_plot_pred, 2, mean)
lnumb_ci<-apply(lnumb_pp_plot_pred, 2, quantile, probs=c(0.05, 0.95))

windows()
par(mfrow=c(2,3)) 

plot(c(0,max(av$ap,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Experienced") 
lines(ap_range, ap_mean, lwd=3, lty=1)
lines(ap_range, ap_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(ap_range, ap_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bp_range, bp_mean, lwd=3, lty=2)
lines(bp_range, bp_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bp_range, bp_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$aally,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Proportion Allied") 
lines(aally_range, aally_mean, lwd=3, lty=1)
lines(aally_range, aally_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(aally_range, aally_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bally_range, bally_mean, lwd=3, lty=2)
lines(bally_range, bally_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bally_range, bally_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$aq,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Previous Collaboration") 
lines(aq_range, aq_mean, lwd=3, lty=1)
lines(aq_range, aq_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(aq_range, aq_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bq_range, bq_mean, lwd=3, lty=2)
lines(bq_range, bq_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bq_range, bq_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$qsuccess,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Collaborative Success Rate") 
lines(qsuccess_range, qsuccess_mean, lwd=3, lty=1)
lines(qsuccess_range, qsuccess_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(qsuccess_range, qsuccess_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bqsuccess_range, bqsuccess_mean, lwd=3, lty=2)
lines(bqsuccess_range, bqsuccess_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bqsuccess_range, bqsuccess_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$ademvs,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Political Variance") 
lines(ademvs_range, ademvs_mean, lwd=3, lty=1)
lines(ademvs_range, ademvs_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(ademvs_range, ademvs_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(bdemvs_range, bdemvs_mean, lwd=3, lty=2)
lines(bdemvs_range, bdemvs_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(bdemvs_range, bdemvs_ci[2,],  lwd=1.5, lty=2, col="grey50")

plot(c(0,max(av$l.numa,na.rm=T)*1.1), c(0,1), xlab="", ylab="Pr (Victory)", type="n", xaxs="i", yaxs="i", main="Log Coalition Size") 
lines(lnuma_range, lnuma_mean, lwd=3, lty=1)
lines(lnuma_range, lnuma_ci[1,], lwd=1.5, lty=1, col="grey50")
lines(lnuma_range, lnuma_ci[2,],  lwd=1.5, lty=1, col="grey50")
lines(lnumb_range, lnumb_mean, lwd=3, lty=2)
lines(lnumb_range, lnumb_ci[1,], lwd=1.5, lty=2, col="grey50")
lines(lnumb_range, lnumb_ci[2,],  lwd=1.5, lty=2, col="grey50")



####### Cleaning up R's memory to make room for Robustness Checks
rm(aally_ci, aally_pp_plot, aally_pp_plot_pred, ademvs_ci, ademvs_pp_plot, ademvs_pp_plot_pred,
            All_sam, All_sam_def, ap_ci, ap_pp_plot, ap_pp_plot_pred, aq_ci, aq_pp_plot, aq_pp_plot_pred, awin_ci, awin_pp_plot, awin_pp_plot_pred,
            bally_ci, bally_pp_plot, bally_pp_plot_pred, bdemvs_ci, bdemvs_pp_plot, bdemvs_pp_plot_pred, bp_ci, bp_pp_plot, bp_pp_plot_pred,
            bq_ci, bq_pp_plot, bq_pp_plot_pred, bqsuccess_ci, bqsuccess_pp_plot, bqsuccess_pp_plot_pred, bwin_ci, bwin_pp_plot, bwin_pp_plot_pred,
            lnuma_ci, lnuma_pp_plot, lnuma_pp_plot_pred, lnumb_ci, lnumb_pp_plot, lnumb_pp_plot_pred, qsuccess_ci, qsuccess_pp_plot, qsuccess_pp_plot_pred)

########## Models discussed in only in footnotes: not logging size of coalition ##########
## Table 1 when size of coalition is not logged
coal_size <- MCMClogit(aout_vic ~ ap + awin  + aally + aq + ademvs + exp(lnuma) + caprat2 + exp(lnumb) + adistm, burnin=1000000, mcmc=5000000, data=coal)
coal_size2 <- MCMClogit(aout_vic ~ ap  + aally + aq + qsuccess +  ademvs + exp(lnuma) + caprat2 + exp(lnumb) + adistm, burnin=1000000, mcmc=5000000, data=coal)

## Table 2 Initiating Coalitions when size of coalition is not logged
asab_size <- MCMClogit(av_out_all ~ ap + awin + aally + aq +  ademvs+ exp(l.numa) + caprat2 + exp(l.numb)+ adistm, burnin=1000000, mcmc=5000000, data=av)
asab_size2 <- MCMClogit(av_out_all ~ ap +  aally + aq + qsuccess + ademvs+ exp(l.numa) + caprat2 + exp(l.numb)+ adistm, burnin=1000000, mcmc=5000000, data=av)

## Table 2 Defending Coalitions when size of coalition is not logged
asdb_size <- MCMClogit(dv_out_all ~ bp + bwin + bally + bq +  bdemvs+ exp(l.numa) + caprat2 + exp(l.numb)+ bdistm, burnin=1000000, mcmc=5000000, data=dv)
asdb_size2 <- MCMClogit(dv_out_all ~ bp +  bally + bq + qsuccess + bdemvs+ exp(l.numa) + caprat2 + exp(l.numb)+ bdistm, burnin=1000000, mcmc=5000000, data=dv)

rm(coal_size, coal_size2, asab_size, asab_size2, asdb_size, asdb_size2)

############################### Appendix Tables & Figures ###############################

##### Model Fit Plots

## ROC and PR Plots Table 1
All_sam<-coal_est[sample(nrow(coal_est),size=10000,replace=F),]
X <- cbind(1, coal[,c("ap", "awin", "aally", "aq", "ademvs", "lnuma", "caprat2","lnumb",  "adistm")])

All_sam2<-coal_est2[sample(nrow(coal_est2),size=10000,replace=F),]
X2 <- cbind(1, coal[,c("ap", "aally", "aq", "qsuccess", "ademvs", "lnuma", "caprat2","lnumb",  "adistm")])

all.pp <- pprobs(All_sam, X)
all.pmeans <- apply(all.pp, 1, mean)
all2.pp <- pprobs(All_sam2, X2)
all2.pmeans <- apply(all2.pp, 1, mean)

###ROC Plot
preds <- prediction(as.numeric(all.pmeans), as.numeric(coal$aout_vic))
perf <- performance(preds, measure="tpr", x.measure="fpr") 
preds2 <- prediction(as.numeric(all2.pmeans), as.numeric(coal$aout_vic))
perf2 <- performance(preds2, measure="tpr", x.measure="fpr") 
performance(preds, "auc")
performance(preds2, "auc")

windows()
plot(perf, main="", lwd=4, col="black")
plot(perf2, add=TRUE, lwd=4, col="black", lty=2)
text(0.6,0.5, "Left Model AUC = 0.809", cex=1.7, col="black")
text(0.6,0.4, "Right Model AUC = 0.795", cex=1.7, col="black")
rug(jitter(coal$aout_vic))

###PR plots 
perf_pr  <- performance(preds, measure="prec", x.measure="rec") #PR
perf2_pr  <- performance(preds2, measure="prec", x.measure="rec") #PR

pr<-pr.curve(scores.class0 = perf_pr@x.values[[1]][2:258], scores.class1 = perf_pr@y.values[[1]][2:258], curve=T)
pr2<-pr.curve(scores.class0 = perf2_pr@x.values[[1]][2:258], scores.class1 = perf2_pr@y.values[[1]][2:258], curve=T)

windows()
plot(pr, main="", lwd=4, col="black",auc.main=FALSE)
plot(pr2, add=T, main="   ", sub="", lwd=4, col="black", lty=2,auc.main=FALSE)
rug(jitter(coal$aout_vic))
text(0.6,0.5, "Left Model AUC = 0.899", cex=1.7, col="black")
text(0.6,0.4, "Right Model AUC = 0.904", cex=1.7, col="black")


## ROC and PR Plots Table 2 Initiating Coalitions
init_sam<-asab_all[sample(nrow(asab_all),size=10000,replace=F),]
X_av <- cbind(1, av[,c("ap", "awin", "aally", "aq", "ademvs", "l.numa", "caprat2","l.numb",  "adistm")])
init_sam2<-asab_all2[sample(nrow(asab_all2),size=10000,replace=F),]
X2_av <- cbind(1, av[,c("ap", "aally", "aq", "qsuccess", "ademvs", "l.numa", "caprat2","l.numb",  "adistm")])

init.pp <- pprobs(init_sam, X_av)
init.pmeans <- apply(init.pp, 1, mean)
init2.pp <- pprobs(init_sam2, X2_av)
init2.pmeans <- apply(init2.pp, 1, mean)

### ROC Plot
preds_init <- prediction(as.numeric(init.pmeans), as.numeric(av$av_out_all))
perf_init <- performance(preds_init, measure="tpr", x.measure="fpr") 
preds_init2 <- prediction(as.numeric(init2.pmeans), as.numeric(av$av_out_all))
perf_init2 <- performance(preds_init2, measure="tpr", x.measure="fpr") 
performance(preds_init, "auc")
performance(preds_init2, "auc")

### Initiating ROC
windows()
plot(perf_init, main="", lwd=4, col="black")
plot(perf_init2, add=TRUE, lwd=4, col="black", lty=2)
text(0.6,0.5, "Left Model AUC = 0.857", cex=1.7, col="black")
text(0.6,0.4, "Right Model AUC = 0.866", cex=1.7, col="black")
rug(jitter(av$av_out_all))

### Initiating PR
perf_init_pr  <- performance(preds_init, measure="prec", x.measure="rec") #PR
perf_init2_pr  <- performance(preds_init2, measure="prec", x.measure="rec") #PR

pr_init<-pr.curve(scores.class0 = perf_init_pr@x.values[[1]][2:127], scores.class1 = perf_pr@y.values[[1]][2:127], curve=T)
pr2_init<-pr.curve(scores.class0 = perf_init2_pr@x.values[[1]][2:127], scores.class1 = perf_init2_pr@y.values[[1]][2:127], curve=T)

windows()
plot(pr_init, main="", lwd=4, col="black",auc.main=FALSE)
plot(pr2_init, add=T, main="   ", sub="", lwd=4, col="black", lty=2,auc.main=FALSE)
rug(jitter(av$av_out_all))
text(0.5,0.5, "Left Model AUC = 0.867", cex=1.7, col="black")
text(0.5,0.4, "Right Model AUC = 0.861", cex=1.7, col="black")

### ROC and PR Plots Table 2 Defending Coalitions
def_sam<-asdb_all[sample(nrow(asdb_all),size=10000,replace=F),]
X_def <- cbind(1, dv[,c("bp", "bwin", "bally", "bq", "bdemvs", "l.numa", "caprat2","l.numb",  "bdistm")])
def_sam2<-asdb_all2[sample(nrow(asdb_all2),size=10000,replace=F),]
X2_def <- cbind(1, dv[,c("bp", "bally", "bq", "qsuccess", "bdemvs", "l.numa", "caprat2","l.numb",  "bdistm")])

def.pp <- pprobs(def_sam, X_def)
def.pmeans <- apply(def.pp, 1, mean)
def2.pp <- pprobs(def_sam2, X2_def)
def2.pmeans <- apply(def2.pp, 1, mean)

### ROC Defending
preds_def <- prediction(as.numeric(def.pmeans), as.numeric(dv$dv_out_all))
perf_def <- performance(preds_def, measure="tpr", x.measure="fpr") 
preds_def2 <- prediction(as.numeric(def2.pmeans), as.numeric(dv$dv_out_all))
perf_def2 <- performance(preds_def2, measure="tpr", x.measure="fpr") 
performance(preds_def, "auc")
performance(preds_def2, "auc")

windows()
plot(perf_def, main="", lwd=4, col="black")
plot(perf_def2, add=TRUE, lwd=4, col="black", lty=2)
text(0.6,0.5, "Left Model AUC = 0.820", cex=1.7, col="black")
text(0.6,0.4, "Right Model AUC = 0.802", cex=1.7, col="black")
rug(jitter(dv$dv_out_all))

### PR Defending
perf_def_pr  <- performance(preds_def, measure="prec", x.measure="rec") #PR
perf_def2_pr  <- performance(preds_def2, measure="prec", x.measure="rec") #PR

pr_def<-pr.curve(scores.class0 = perf_def_pr@x.values[[1]][2:169], scores.class1 = perf_pr@y.values[[1]][2:169], curve=T)
pr2_def<-pr.curve(scores.class0 = perf_def2_pr@x.values[[1]][2:169], scores.class1 = perf_def2_pr@y.values[[1]][2:169], curve=T)

windows()
plot(pr_def, main="", lwd=4, col="black",auc.main=FALSE)
plot(pr2_def, add=T, main="   ", sub="", lwd=4, col="black", lty=2,auc.main=FALSE)
rug(jitter(dv$dv_out_all))
text(0.5,0.5, "Left Model AUC = 0.894", cex=1.7, col="black")
text(0.5,0.4, "Right Model AUC = 0.887", cex=1.7, col="black")

rug(jitter(dv$dv_out_all))

### Separation Plots
## Table 1
windows()
par(mfrow = c(2, 1))
separationplot(all.pmeans, coal$aout_vic, 
               heading = "Left Model in Table 1", newplot = FALSE)
separationplot(all2.pmeans, coal$aout_vic, 
               heading = "Right Model in Table 1", newplot = FALSE)

## Table 2
windows()
par(mfrow = c(4, 1))
separationplot(init.pmeans, av$av_out_all, 
               heading = "Left Initiating Coal. Model in Table 2", newplot = FALSE)
separationplot(init2.pmeans, av$av_out_all, 
               heading = "Right Initiating Coal. Model in Table 2", newplot = FALSE)
separationplot(def.pmeans, dv$dv_out_all, 
               heading = "Left Defending Coal. Model in Table 2", newplot = FALSE)
separationplot(def2.pmeans, dv$dv_out_all, 
               heading = "Right Defending Coal. Model in Table 2", newplot = FALSE)


rm(All_sam, All_sam2, all.pp, all2.pp, def_sam, def_sam2, def.pp, def2.pp, init_sam, init_sam2, init.pp, init2.pp, X, X2, X_av, X_def, X2_av, X2_def)
rm(perf, perf_def, perf_def_pr, perf_def2, perf_def2_pr, perf_init, perf_init_pr, perf_init2, perf_init2_pr, 
            perf_pr, perf2, perf2_pr, pr, pr_def, pr_init, pr2, pr2_def, pr2_init, preds, preds_def, preds_def2,
            preds_init, preds_init2, preds2)

###### Descriptive Statistics
xtable(rbind(summary(coal$aout_vic),summary(coal$ap),summary(coal$awin),summary(coal$aally),summary(coal$aq),summary(coal$qsuccess),summary(coal$ademvs),summary(coal$lnuma),summary(coal$caprat2),summary(coal$lnumb),summary(coal$adistm)))
xtable(rbind(summary(av$av_out_all),summary(av$ap),summary(av$awin),summary(av$aally),summary(av$aq),summary(av$qsuccess),summary(av$ademvs),summary(av$l.numa),summary(av$caprat2),summary(av$l.numb),summary(av$adistm)))
xtable(rbind(summary(dv$dv_out_all),summary(dv$bp),summary(dv$bwin),summary(dv$bally),summary(dv$bq),summary(dv$qsuccess),summary(dv$bdemvs),summary(dv$l.numa),summary(dv$caprat2),summary(dv$l.numb),summary(dv$bdistm)))

###### Alternative specifications of Table 1 in Appendix:
## With both Success Rate and Collaborative Success Rate in the same model
coal_est_app1 <- MCMClogit(aout_vic ~ ap + awin  + aally + aq + qsuccess +  ademvs + lnuma + caprat2 + lnumb + adistm, burnin=1000000, mcmc=5000000, data=coal)
xtable(interpret(coal_est_app1))

## Controlling for dispute issue
coal_est_app2 <- MCMClogit(aout_vic ~ ap + awin  + aally + aq +  ademvs + lnuma + caprat2 + lnumb + adistm + ter + pol, burnin=1000000, mcmc=5000000, data=coal)
coal_est_app2b <- MCMClogit(aout_vic ~ ap +  aally + aq +  qsuccess + ademvs + lnuma + caprat2 + lnumb + adistm + ter + pol, burnin=1000000, mcmc=5000000, data=coal)

xtable(cbind(tex.interpret(coal_est_app2),tex.interpret(coal_est_app2b)))

### Coal v single state
coalvsing<-subset(coal, lnumb==0)

coal_sing <- MCMClogit(aout_vic ~ ap + awin + aally + aq + ademvs + lnuma + caprat2 + adistm, burnin=1000000, mcmc=10000000, data=coalvsing)
coal_sing2 <- MCMClogit(aout_vic ~ ap + aally + aq + qsuccess+ademvs + lnuma + caprat2 + adistm, burnin=1000000, mcmc=10000000, data=coalvsing)
xtable(cbind(tex.interpret(coal_sing),tex.interpret(coal_sing2)))

### Coalition v coalition
coalvcoal<-subset(coal, lnumb>0)

## Variables capture difference between coalitions
diffp<-coalvcoal$ap-coalvcoal$bp
diffwin<-coalvcoal$awin-coalvcoal$bwin
diffdemvs<-coalvcoal$ademvs-coalvcoal$bdemvs
diffally<-coalvcoal$aally-coalvcoal$bally
diffq<-coalvcoal$aq-coalvcoal$bq
difflnum<-coalvcoal$lnuma-coalvcoal$lnumb
diffdistm<-coalvcoal$adistm-coalvcoal$bdistm
diffqsuccess<-coalvcoal$qsuccess-coalvcoal$qsuccessb
coalvcoal<-cbind(coalvcoal, diffp, diffwin, diffdemvs, diffally, diffq, difflnum, diffdistm, diffqsuccess)

coal_coal <- MCMCoprobit(aoutt~diffp + diffwin +  diffally + diffq+ diffdemvs+difflnum + caprat2 + diffdistm, burnin=1000000, mcmc=5000000, data=coalvcoal)
coal_coal2 <- MCMCoprobit(aoutt~diffp +  diffally + diffq+ diffqsuccess+diffdemvs+difflnum + caprat2 + diffdistm, burnin=1000000, mcmc=5000000, data=coalvcoal)

xtable(cbind(interpret(coal_coal), interpret(coal_coal2)))

rm(coalvsing, coalvcoal, coal_est_app1, coal_est_app2, coal_est_app2b, coal_sing, coal_sing2, coal_coal, coal_coal2)
 
#### Table 1 with All MIDS
all<-read.csv("Replication Files/allMIDs.csv")

coal_all <- MCMClogit(formcoal, burnin=1000000, mcmc=5000000, data=all)
coal_all2 <- MCMClogit(aout_vic ~ ap  + aally + aq + qsuccess_all +  ademvs + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=all)

xtable(cbind(interpret(coal_all), interpret(coal_all2)))

### MLE Table 1
coal_mle <- glm(formcoal, data=coal, family=binomial(link="logit"))
coal_mle2 <- glm(formcoal2, data=coal, family=binomial(link="logit"))
xtable(cbind(summary(coal_mle)$coef, summary(coal_mle2)$coef))

rm(coal_all, coal_all2, coal_mle, coal_mle2)



###### Alternative specifications of Table 2 in Appendix:
## With both Success Rate and Collaborative Success Rate in the same model
asab_est3 <- MCMClogit(av_out_all ~ ap + awin + aally + aq + qsuccess + ademvs+ l.numa + caprat2 + l.numb+ adistm, burnin=1000000, mcmc=5000000, data=av)
asdb_est3 <- MCMClogit(dv_out_all ~ bp + bwin + bally + bq +  qsuccess+bdemvs+ l.numa + caprat2 + l.numb+ bdistm, burnin=1000000, mcmc=5000000, data=dv)

xtable(cbind(tex.interpret(asab_est3), tex.interpret(asdb_est3)))

## Controlling for dispute issue
av_app1 <- MCMClogit(av_out_all ~ ap + awin + aally + aq +  ademvs+ l.numa + caprat2 + l.numb+ adistm+ter+pol, burnin=1000000, mcmc=5000000, data=av)
av_app2 <- MCMClogit(av_out_all ~ ap +  aally + aq + qsuccess+ ademvs+ l.numa + caprat2 + l.numb+ adistm+ter+pol, burnin=1000000, mcmc=5000000, data=av)

dv_app1 <- MCMClogit(dv_out_all ~ bp + bwin + bally + bq +  bdemvs+ l.numa + caprat2 + l.numb+ bdistm+ter+pol, burnin=1000000, mcmc=5000000, data=dv)
dv_app2 <- MCMClogit(dv_out_all ~ bp + bally + bq +  qsuccess+ bdemvs+ l.numa + caprat2 + l.numb+ bdistm+ter+pol, burnin=1000000, mcmc=5000000, data=dv)

xtable(cbind(tex.interpret(av_app1), tex.interpret(av_app2)))
xtable(cbind(tex.interpret(dv_app1), tex.interpret(dv_app2)))

rm(asab_est3, asdb_est3, av_app1, av_app2, dv_app1, dv_app2)

#### Table 2 with All MIDS
## Initiating
av_all <- subset(subset(all, numa > 1))

coal_av_all <- MCMClogit(aout_vic ~ ap + awin  + aally + aq +   ademvs + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=av_all)
coal_av_all2 <- MCMClogit(aout_vic ~ ap  + aally + aq + qsuccessa_all +  ademvs + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=av_all)

xtable(cbind(interpret(coal_av_all), interpret(coal_av_all2)))

## Defending
dv_all <- subset(subset(all, numb > 1))
dv_all$bout_vic<-ifelse(dv_all$aout_vic==1, 0, 1)

coal_dv_all <- MCMClogit(bout_vic ~ bp + bwin  + bally + bq +   bdemvs + lnuma + caprat2 + lnumb+ bdistm, burnin=1000000, mcmc=5000000, data=dv_all)
coal_dv_all2 <- MCMClogit(bout_vic ~ bp  + bally + bq + qsuccessb_all +  bdemvs + lnuma + caprat2 + lnumb+ bdistm, burnin=1000000, mcmc=5000000, data=dv_all)

xtable(cbind(interpret(coal_dv_all), interpret(coal_dv_all2)))


####MLE Estimates
asab_mle <- glm(formBaall1, family=binomial(link="logit"), data=av)
asab_mle2 <- glm(formBaall2, family=binomial(link="logit"), data=av)
xtable(cbind(summary(asab_mle)$coef, summary(asab_mle2)$coef))

asdb_mle <- glm(formBdall1, family=binomial(link="logit"), data=dv)
asdb_mle2 <- glm(formBdall2, family=binomial(link="logit"), data=dv)
xtable(cbind(summary(asdb_mle)$coef, summary(asdb_mle2)$coef))

rm(av_all, dv_all, coal_av_all, coal_av_all2, coal_dv_all, coal_dv_all2, asab_mle, asab_mle2, asdb_mle, asdb_mle2)

####### Alternative measures of legitimacy
#### UN Ideal Points
## Ignoring observations with varid missing; UN Ideal Points only available in part of the time series
coal_id<-subset(coal, varid>=0)


All_IDpt <- MCMClogit(aout_vic ~ ap + awin  + aally + aq + varid + lnuma + caprat2 + lnumb+adistm, burnin=1000000, mcmc=5000000, data=coal_id)
All_IDpt2 <- MCMClogit(aout_vic ~ ap + aally + aq+ qsuccess + varid +  lnuma + caprat2 +lnumb+adistm, burnin=1000000, mcmc=5000000, data=coal_id)

xtable(cbind(tex.interpret(All_IDpt), tex.interpret(All_IDpt2)))


#### Global & Regional support
coal_UNsupport <- MCMClogit(aout_vic ~ ap +awin + aally+ aq+ unsupportstrong + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=coal)
coal_UNsupport2 <- MCMClogit(aout_vic ~ ap + aally+ aq+ qsuccess + unsupportstrong + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=coal)

xtable(cbind(interpret(coal_UNsupport), interpret(coal_UNsupport2)))

unregstrong<-vector(length=nrow(coal))
for(i in 1:nrow(coal)){
  if(is.na(coal$unsupportstrong[i])=="TRUE"){unregstrong[i]<-"NA"}
  if(is.na(coal$unsupportstrong[i])!="TRUE"){
    if(coal$unsupportstrong[i]==1 | coal$regsupportstrong[i]==1){unregstrong[i]<-1}
    if(coal$unsupportstrong[i]==0 & coal$regsupportstrong[i]==0){unregstrong[i]<-0}
  }
}
unregstrong<-as.numeric(unregstrong)
coal<-cbind(coal, unregstrong)


coal_UN_regsupport <- MCMClogit(aout_vic ~ ap +awin + aally+ aq+ unregstrong + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=coal)
coal_UN_regsupport2 <- MCMClogit(aout_vic ~ ap + aally+ aq+ qsuccess + unregstrong + lnuma + caprat2 + lnumb+ adistm, burnin=1000000, mcmc=5000000, data=coal)

xtable(cbind(interpret(coal_UN_regsupport), interpret(coal_UN_regsupport2)))

rm(coal_id, All_IDpt, All_IDpt2, coal_UN_regsupport2, coal_UN_regsupport, coal_UNsupport2, coal_UNsupport)


###### Additional robustness checks
### Robustness Check: Outliers
coalr <- subset(coal, dispnum != c(4137))

coal_out <- MCMClogit(formcoal, burnin=1000000, mcmc=5000000, data=coalr)
coal_out2 <- MCMClogit(formcoal2, burnin=1000000, mcmc=5000000, data=coalr)

xtable(cbind(interpret(coal_out), interpret(coal_out2)))


### Initiating and Defending
avr <- subset(av, dispnum != c(4137))
dvr <- subset(dv, dispnum != c(3957))

asab_out <- MCMClogit(formBaall1, burnin=1000000, mcmc=5000000, data=avr)
asab2_out <- MCMClogit(formBaall2, burnin=1000000, mcmc=5000000, data=avr)
xtable(cbind(interpret(asab_out), interpret(asab2_out)))

asdb_out <- MCMClogit(formBdall1, burnin=1000000, mcmc=5000000, data=dvr)
asdb2_out <- MCMClogit(formBdall2, burnin=1000000, mcmc=5000000, data=dvr)
xtable(cbind(interpret(asdb_out), interpret(asdb2_out)))

rm(coal_out, coal_out2, asab_out, asab2_out, asdb_out, asdb2_out)

### Robustness Check: Temporal Heterogeneity
# Start Year
coal_st <- MCMClogit(aout_vic ~ ap +awin + aally+ aq+ ademvs + lnuma + caprat2 + lnumb+ adistm+styear, burnin=1000000, mcmc=5000000, data=coal)
coal_st2 <- MCMClogit(aout_vic ~ ap + aally+ aq+ qsuccess + ademvs + lnuma + caprat2 + lnumb+ adistm+styear, burnin=1000000, mcmc=5000000, data=coal)

xtable(cbind(interpret(coal_st), interpret(coal_st2)))

# Post WWII

postwwii<-vector()
for(i in 1:nrow(coal)){
  if(coal$styear[i]<=1945){postwwii[i]<-0}
  if(coal$styear[i]>1945){postwwii[i]<-1}
}

coal_post <- MCMClogit(aout_vic ~ ap +awin + aally+ aq+ ademvs + lnuma + caprat2 + lnumb+ adistm+postwwii, burnin=1000000, mcmc=5000000, data=coal)
coal_post2 <- MCMClogit(aout_vic ~ ap + aally+ aq+ qsuccess + ademvs + lnuma + caprat2 + lnumb+ adistm+postwwii, burnin=1000000, mcmc=5000000, data=coal)

xtable(cbind(interpret(coal_post), interpret(coal_post2)))

rm(coal_st, coal_st2, coal_post, coal_post2)

### Initiating and Defending
# Start Year
coal_av_st <- MCMClogit(av_out_all ~ ap + awin + aally + aq +  ademvs+ l.numa + caprat2 + l.numb+ adistm+styear, burnin=1000000, mcmc=5000000, data=av)
coal_av_st2 <- MCMClogit(av_out_all ~ ap +  aally + aq + qsuccess + ademvs+ l.numa + caprat2 + l.numb+ adistm+styear, burnin=1000000, mcmc=5000000, data=av)

xtable(cbind(interpret(coal_av_st), interpret(coal_av_st2)))

coal_dv_st <- MCMClogit(dv_out_all ~ bp + bwin + bally + bq +  bdemvs+ l.numa + caprat2 + l.numb+ bdistm+styear, burnin=1000000, mcmc=5000000, data=dv)
coal_dv_st2 <- MCMClogit(dv_out_all ~ bp +  bally + bq + qsuccess + bdemvs+ l.numa + caprat2 + l.numb+ bdistm+styear, burnin=1000000, mcmc=5000000, data=dv)

xtable(cbind(interpret(coal_dv_st), interpret(coal_dv_st2)))

# Post WWII

apostwwii<-vector()
for(i in 1:nrow(av)){
  if(av$styear[i]<=1945){apostwwii[i]<-0}
  if(av$styear[i]>1945){apostwwii[i]<-1}
}

dpostwwii<-vector()
for(i in 1:nrow(dv)){
  if(dv$styear[i]<=1945){dpostwwii[i]<-0}
  if(dv$styear[i]>1945){dpostwwii[i]<-1}
}

coal_av_post <- MCMClogit(av_out_all ~ ap + awin + aally + aq +  ademvs+ l.numa + caprat2 + l.numb+ adistm+apostwwii, burnin=1000000, mcmc=5000000, data=av)
coal_av_post2 <- MCMClogit(av_out_all ~ ap +  aally + aq + qsuccess + ademvs+ l.numa + caprat2 + l.numb+ adistm+apostwwii, burnin=1000000, mcmc=5000000, data=av)

xtable(cbind(interpret(coal_av_post), interpret(coal_av_post2)))

coal_dv_post <-MCMClogit(dv_out_all ~ bp + bwin + bally + bq +  bdemvs+ l.numa + caprat2 + l.numb+ bdistm+dpostwwii, burnin=1000000, mcmc=5000000, data=dv)
coal_dv_post2 <- MCMClogit(dv_out_all ~ bp +  bally + bq + qsuccess + bdemvs+ l.numa + caprat2 + l.numb+ bdistm+dpostwwii, burnin=1000000, mcmc=5000000, data=dv)

xtable(cbind(interpret(coal_dv_post), interpret(coal_dv_post2)))

rm(coal_av_post, coal_av_post2, coal_dv_post, coal_dv_post2, coal_av_st, coal_av_st2, coal_dv_st, coal_dv_st2)

rm(all, avr, coalr, dvr, aally_mean, aally_range, ademvs_mean, ademvs_range, all.pmeans, all2.pmeans, ap_mean, ap_range, apostwwii, aq_mean,
            aq_range, awin_mean, awin_range, bally_mean, bally_range, bdemvs_mean, bdemvs_range, bp_mean, bp_range, bq_mean, bq_range,
            bqsuccess_mean, bqsuccess_range, bwin_mean, bwin_range, def.pmeans, def2.pmeans, diffally, diffdemvs, diffdistm, difflnum, diffp, 
            diffq, diffqsuccess, diffwin, dpostwwii, i, init.pmeans, init2.pmeans, lnuma_mean, lnuma_range, lnumb_mean, lnumb_range,
            postwwii, qsuccess_mean, qsuccess_range, unregstrong, pprobs)

#### Convergence Diagnostics for Tables in Main Text
### Table 1
windows()
par(mfrow=c(2, 5))
traceplot(coal_est[,2], main="Proportion Experienced")
traceplot(coal_est[,3], main="Success Rate")
traceplot(coal_est[,4], main="Proportion Allied")
traceplot(coal_est[,5], main="Previous Collaboration")
traceplot(coal_est[,6], main="Polity Variance")
traceplot(coal_est[,7], main="Log Coalition Size")
traceplot(coal_est[,8], main="Capability Ratio")
traceplot(coal_est[,9], main="Size Side B")
traceplot(coal_est[,10], main="Mean Geog. Distance")
traceplot(coal_est[,1], main="Intercept")

windows()
par(mfrow=c(2, 5))
traceplot(coal_est2[,2], main="Proportion Experienced")
traceplot(coal_est2[,3], main="Proportion Allied")
traceplot(coal_est2[,4], main="Previous Collaboration")
traceplot(coal_est2[,5], main="Collaborative Success Rate")
traceplot(coal_est2[,6], main="Polity Variance")
traceplot(coal_est2[,7], main="Log Coalition Size")
traceplot(coal_est2[,8], main="Capability Ratio")
traceplot(coal_est2[,9], main="Size Side B")
traceplot(coal_est2[,10], main="Mean Geog. Distance")
traceplot(coal_est2[,1], main="Intercept")

windows()
par(mfrow=c(2,5))
autocorr.plot(coal_est[,2], main="Proportion Experienced", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,3], main="Success Rate", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,4], main="Proportion Allied", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,5], main="Previous Collaboration", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,6], main="Polity Variance", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,7], main="Log Coalition Size", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,8], main="Capability Ratio", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,9], main="Size Side B", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,10], main="Mean Geog. Distance", lag.max=100, auto.layout=F)
autocorr.plot(coal_est[,1], main="Intercept", lag.max=100, auto.layout=F)

windows()
par(mfrow=c(2, 5))
autocorr.plot(coal_est2[,2], main="Proportion Experienced", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,3], main="Proportion Allied", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,4], main="Previous Collaboration", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,5], main="Collaborative Success Rate", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,6], main="Polity Variance", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,7], main="Log Coalition Size", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,8], main="Capability Ratio", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,9], main="Size Side B", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,10], main="Mean Geog. Distance", lag.max=100, auto.layout=F)
autocorr.plot(coal_est2[,1], main="Intercept", lag.max=100, auto.layout=F)


ab<-geweke.diag(coal_est, frac1=.25)
db<-geweke.diag(coal_est2, frac1=.25)
xtable(as.matrix(cbind(ab$z, db$z)))

rm(ab, db)


#Gelman Rubin
coal_est3 <- MCMClogit(formcoal, burnin=1000000, mcmc=10000000, data=coal, beta.start=-3)
coal_est4 <- MCMClogit(formcoal, burnin=1000000, mcmc=10000000, data=coal, beta.start=3)

gelman.diag(as.mcmc.list(list(coal_est,  coal_est3, coal_est4)), autoburnin = FALSE)

rm(coal_est3, coal_est4)

coal_est5 <- MCMClogit(formcoal2, burnin=1000000, mcmc=10000000, data=coal, beta.start=-3)
coal_est6 <- MCMClogit(formcoal2, burnin=1000000, mcmc=10000000, data=coal, beta.start=3)

gelman.diag(as.mcmc.list(list(coal_est2, coal_est5, coal_est6)), autoburnin = FALSE)

rm(coal_est5, coal_est6)

#Heidelberger & Welch
heidel.diag(coal_est)
heidel.diag(coal_est2)

rm(coal, coal_est, coal_est2)

##Diagnostics for Table 2
windows()
par(mfrow=c(2, 5))
traceplot(asab_all[,2], main="Proportion Experienced")
traceplot(asab_all[,3], main="Success Rate")
traceplot(asab_all[,4], main="Proportion Allied")
traceplot(asab_all[,5], main="Previous Collaboration")
traceplot(asab_all[,6], main="Polity Variance")
traceplot(asab_all[,7], main="Log Coalition Size")
traceplot(asab_all[,8], main="Capability Ratio")
traceplot(asab_all[,9], main="Size Side B")
traceplot(asab_all[,10], main="Mean Geog. Distance")
traceplot(asab_all[,1], main="Intercept")

windows()
par(mfrow=c(2, 5))
traceplot(asab_all2[,2], main="Proportion Experienced")
traceplot(asab_all2[,3], main="Proportion Allied")
traceplot(asab_all2[,4], main="Previous Collaboration")
traceplot(asab_all2[,5], main="Collaborative Success Rate")
traceplot(asab_all2[,6], main="Polity Variance")
traceplot(asab_all2[,7], main="Log Coalition Size")
traceplot(asab_all2[,8], main="Capability Ratio")
traceplot(asab_all2[,9], main="Size Side B")
traceplot(asab_all2[,10], main="Mean Geog. Distance")
traceplot(asab_all2[,1], main="Intercept")

windows()
par(mfrow=c(2, 5))
traceplot(asdb_all[,2], main="Proportion Experienced")
traceplot(asdb_all[,3], main="Success Rate")
traceplot(asdb_all[,4], main="Proportion Allied")
traceplot(asdb_all[,5], main="Previous Collaboration")
traceplot(asdb_all[,6], main="Polity Variance")
traceplot(asdb_all[,7], main="Log Coalition Size")
traceplot(asdb_all[,8], main="Capability Ratio")
traceplot(asdb_all[,9], main="Size Side B")
traceplot(asdb_all[,10], main="Mean Geog. Distance")
traceplot(asdb_all[,1], main="Intercept")

windows()
par(mfrow=c(2, 5))
traceplot(asdb_all2[,2], main="Proportion Experienced")
traceplot(asdb_all2[,3], main="Proportion Allied")
traceplot(asdb_all2[,4], main="Previous Collaboration")
traceplot(asdb_all2[,5], main="Collaborative Success Rate")
traceplot(asdb_all2[,6], main="Polity Variance")
traceplot(asdb_all2[,7], main="Log Coalition Size")
traceplot(asdb_all2[,8], main="Capability Ratio")
traceplot(asdb_all2[,9], main="Size Side B")
traceplot(asdb_all2[,10], main="Mean Geog. Distance")
traceplot(asdb_all2[,1], main="Intercept")

windows()
par(mfrow=c(2,5))
autocorr.plot(asab_all[,2], main="Proportion Experienced", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,3], main="Success Rate", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,4], main="Proportion Allied", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,5], main="Previous Collaboration", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,6], main="Polity Variance", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,7], main="Log Coalition Size", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,8], main="Capability Ratio", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,9], main="Size Side B", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,10], main="Mean Geog. Distance", lag.max=100, auto.layout=F)
autocorr.plot(asab_all[,1], main="Intercept", lag.max=100, auto.layout=F)

windows()
par(mfrow=c(2, 5))
autocorr.plot(asab_all2[,2], main="Proportion Experienced", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,3], main="Proportion Allied", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,4], main="Previous Collaboration", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,5], main="Collaborative Success Rate", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,6], main="Polity Variance", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,7], main="Log Coalition Size", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,8], main="Capability Ratio", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,9], main="Size Side B", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,10], main="Mean Geog. Distance", lag.max=100, auto.layout=F)
autocorr.plot(asab_all2[,1], main="Intercept", lag.max=100, auto.layout=F)

windows()
par(mfrow=c(2,5))
autocorr.plot(asdb_all[,2], main="Proportion Experienced", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,3], main="Success Rate", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,4], main="Proportion Allied", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,5], main="Previous Collaboration", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,6], main="Polity Variance", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,7], main="Log Coalition Size", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,8], main="Capability Ratio", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,9], main="Size Side B", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,10], main="Mean Geog. Distance", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all[,1], main="Intercept", lag.max=100, auto.layout=F)

windows()
par(mfrow=c(2, 5))
autocorr.plot(asdb_all2[,2], main="Proportion Experienced", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,3], main="Proportion Allied", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,4], main="Previous Collaboration", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,5], main="Collaborative Success Rate", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,6], main="Polity Variance", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,7], main="Log Coalition Size", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,8], main="Capability Ratio", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,9], main="Size Side B", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,10], main="Mean Geog. Distance", lag.max=100, auto.layout=F)
autocorr.plot(asdb_all2[,1], main="Intercept", lag.max=100, auto.layout=F)


a<-geweke.diag(asab_all, frac1=.25)
b<-geweke.diag(asab_all2, frac1=.25)
c<-geweke.diag(asdb_all, frac1=.25)
d<-geweke.diag(asdb_all2, frac1=.25)

xtable(as.matrix(cbind(a$z, b$z, c$z, d$z)))

rm(a, b, c, d)

#Gelman Rubin
asab_all3 <- MCMClogit(formBaall1, burnin=1000000, mcmc=5000000, data=av, beta.start=-3)
asab_all4 <- MCMClogit(formBaall1, burnin=1000000, mcmc=5000000, data=av, beta.start=3)

gelman.diag(as.mcmc.list(list(asab_all,  asab_all3, asab_all4)), autoburnin = FALSE)

rm(asab_all3, asab_all4)

asab_all5 <- MCMClogit(formBaall2, burnin=1000000, mcmc=5000000, data=av, beta.start=-3)
asab_all6 <- MCMClogit(formBaall2, burnin=1000000, mcmc=5000000, data=av, beta.start=3)

gelman.diag(as.mcmc.list(list(asab_all2,  asab_all5, asab_all6)), autoburnin = FALSE)

rm(asab_all5, asab_all6)

####Gelman Rubin for defending coalitions
asdb_all3 <- MCMClogit(formBdall1, burnin=1000000, mcmc=5000000, data=dv, beta.start=-3)
asdb_all4 <- MCMClogit(formBdall1, burnin=1000000, mcmc=5000000, data=dv, beta.start=3)

gelman.diag(as.mcmc.list(list(asdb_all,  asdb_all3, asdb_all4)), autoburnin = FALSE)

rm(asdb_all3, asdb_all4)

asdb_all5 <- MCMClogit(formBdall2, burnin=1000000, mcmc=5000000, data=dv, beta.start=-3)
asdb_all6 <- MCMClogit(formBdall2, burnin=1000000, mcmc=5000000, data=dv, beta.start=3)

gelman.diag(as.mcmc.list(list(asdb_all2,  asdb_all5, asdb_all6)), autoburnin = FALSE)

rm(asdb_all5, asdb_all6)

#Heidelberger & Welch
heidel.diag(asab_all)
heidel.diag(asab_all2)
heidel.diag(asdb_all)
heidel.diag(asdb_all2)